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ABSTRACT 



Aims. We present an elegant method of determining the eigensolutions of the induction and the dynamo equation in a fluid embedded 



in a vacuum. 



Methods. The magnetic field is expanded in a complete set of functions. The new method is based on the biorthogonality of the adjoint 
electric current and the vector potential with an inner product defined by a volume integral over the fluid domain. The advantage of 
this method is that the velocity and the dynamo coefficients of the induction and the dynamo equation do not have to be differentiated 
and thus even numerically determined tabulated values of the coefficients produce reasonable results. 

Results. We provide test calculations and compare with published results obtained by the classical treatment based on the biorthog- 
onality of the magnetic field and its adjoint. We especially consider dynamos with mean-field coefficients determined from direct 
numerical simulations of the geodynamo and compare with initial value calculations and the full MHD simulations. 

Key words. Magnetohydrodynamics (MHD) - magnetic fields - methods: numerical 



1. Introduction 

The generation and evolution of magnetic fields in cosmic bod- 
ies like the planets and stars is generally thought to be governed 
by induction processes due to motions in their electrically con- 
ducting fluid interior. The magnetic field B is described by the 
induction equation 



dB 



at 

where 



= VxDB 



. DB = uxB -iflxB 



(1) 



(2) 



Here u represents the velocity and 77 the magnetic diffusivity. 

In the framework of mean-field theory (e.g. Moffatt 1978; 
Krause & Radler 1980), u and B are considered as mean, e.g. 
ensemble averaged, quantities, whereas the action of the small- 
scale turbulent flow on the mean magnetic field is parametrised 
by the so-called dynamo coefficients, or and 6. They are, in gen- 
eral, tensors of second and third rank, respectively. We use the 
following compact notation of the mean field coefficients, which 
include the so-called y, 6, and K-effects, see e.g. Radler (1980). 
Then, the operator D reads 



DB = uxB + a B - B : (VB)-j]VxB 



(3) 



instead of (2), and acts on the mean magnetic field. Except for 
the additional a and B terms in the D operator, the induction 
and the dynamo equation are formally equivalent. Thus, the new 
method presented here equally applies to both. 

The dynamo region is located in a flow domain V with exte- 
rior vacuum E. In this work we assume V to be either a sphere 



or a spherical shell. The magnetic field B is continuous through 
the boundary dV and potential in E. 

In kinematic dynamo theory, all coefficients (u, a, B, and 77) 
are assumed given and independent of the magnetic field. Thus 
the dynamo equation is linear in the magnetic field and can be 
solved by considering an eigenvalue problem 



AB = VxDB 



(4) 
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with eigenvalues A describing the time evolution proportional to 
exp(At) of the magnetic field B. 

Many studies have been made of the eigenvalues of the 
dynamo operator for various celestial bodies and with many 
forms of the dynamo coefficients (e.g. Bullard & Gellman 1954; 
Roberts 1960; Steenbeck & Krause 1969; Deinzer & Stix 1971; 
Roberts 1972; Roberts & Stix 1972; Gubbins 1973; Kumar & 
Roberts 1975; Schmitt & Schussler 1989; Dudley & James 1989; 
Deinzer et al. 1993; Gubbins et al. 2000; Schubert & Zhang 
2001; Livermore & Jackson 2004, 2005; Jiang & Wang 2006, 
2007). Often the coefficients are approximated by simple analyt- 
ical functions of position, and their tensorial character is disre- 
garded. Recently, the test-field method, developed by Schrinner 
et al. (2005, 2007) (see also Ossendrijver et al. 2001, 2002), al- 
lows one to determine all tensorial components of a and B di- 
rectly from self-consistent numerical simulations (Brandenburg 
et al. 2008; Kapyla et al. 2009). These coefficients are sometimes 
strongly varying functions of position. This may introduce large 
errors because the dynamo operator VxD involves differentia- 
tion of the dynamo coefficients, and these are only available as 
numerically determined tabulated values. 

In this paper we present a new method that does not require 
differentiation, so it is also applicable to numerically determined 
dynamo coefficients. The method is based on the biorthogonality 
of the electric current and the vector potential with an inner prod- 
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uct defined by a volume integral over the fluid. This property has 
already been noted by Radler & Brauer (1987), Hoyng (1988), 
Fuchs et al. (1993), Hoyng & Schutgens (1995), and Radler et 
al. (2002). 

The method is described in detail in Sect. 2. Extensive test 
calculations have been performed and compared with published 
results by other eigenvalue methods. Some of these tests are pre- 
sented in Sect. 3. In Sect. 4 we apply the new method to eigen- 
modes of the dynamo operator with coefficients obtained from 
geodynamo models. Our conclusion are drawn in Sect. 5. 



2. Eigenvalue problem 

We expand the field B of the dynamo in a complete set of func- 
tions A;(r): 

B = Y J e i b i {r) = e i b i . (5) 

i 

Here and in the following we make use of the summation con- 
vention for two identical indices. The expansion functions are 
often eigenfunctions of some differential operator. Since this op- 
erator is, in general, not self-adjoint, the functions are not or- 
thogonal. This problem is handled by using the adjoint set b k (r), 
with the following inner product 

(b k ,bdx= f b k -bid 3 r = 6 ki . (6) 
Jx 

The integration volume X can be either the whole space V + E 
or the fluid domain V alone. The base functions b k (r) and ft,(r) 
constitute a biorthogonal set. For a given set of functions ft,, the 
adjoint set b k depends on the choice of the integration domain 
X, so in principle we have two different sets b k , one for X = V 
and one for X = V + E. 

Later we adopt the free magnetic decay modes, for which 
the base functions ft, and their adjoints b k are known. But at this 
point there is no need to specify which set ft, we actually use. 

2. 1 . Biorthogonal sets 

Starting from a set b k and ft, that is biorthogonal on V + E, a 
very useful biorthogonal set on V is provided by the associated 
electric current j = Vxft and the vector potential a where Vxa = 
ft, with the inner product 

(b k ,bi) v +E = (jk,adv = (a k ,ji)v = hi ■ (7) 

Here we have absorbed a factor of 4n/c in the definition of the 
current j. The relation (7) is derived with the help of the vector 
identity V • (a, x b k ) = b k ■ ( V x a,) - a, • • ( V x b k ) = b k ■ ft, ■■ - a, ■■ ■ j k 
and a volume integration over V + E; a and ft go fast enough 
to zero at infinity. Surface integrals vanish because the field and 
the vector potential are continuous through dV. Volume integrals 
containing currents are restricted to V since j = in E. The inner 
product (7) is invariant under a gauge transformation a — > a+ 
because j^j ■ V^d 3 r = f y V ■ (iffj)d 3 r = f gy i/fj ■ dV = 0, as 
currents and their adjoints run parallel to the boundary. 

Electric currents and vector potentials thus form a biorthog- 
onal set on V. This is essential for the new eigenvalue method 
presented in Sect. 2.3. 



2.2. Classical eigenvalue method 

Inserting the expansion (5) in the dynamo eigenvalue equation 
(4) yields 

Aefii = Vx(De,ft,) . (8) 

Subsequently, we take the inner product (6) based on V with the 
adjoint magnetic field. This leads to 

Ae k = M ki ei with M ki = (b k ,VxDbi) v . (9) 

A partial integration to shift the curl from the second to the first 
term, as done in (7) above and used in the new method below, 
is not possible because the surface term J gv (DbiXb k ) ■ & 2 <t need 
not vanish here. 

We mention as an aside that the magnetic field is often de- 
composed in its poloidal and toroidal components (see Appendix 
A) after which the dynamo equation is formulated in terms of 
the defining scalars P and T. If the dynamo coefficients possess 
certain symmetry properties, the solutions can be split into two 
independent subsets, describing magnetic fields symmetric and 
antisymmetric with respect to the equator. 

2.3. New eigenvalue method 

We start again with (8), which we uncurl to obtain 

Aem, ■ = Dejb, ■ + Vifr . (10) 

Taking now the inner product (7) with the adjoint current results 
in 

Ae k = N ki e t with = (j k ,Dbd v . (11) 

The gradient term drops out as discussed in Sect. 2. 1 above. The 
corresponding adjoint functions b k here are different from those 
in Sect. 2.2 as they pertain to a different inner product. 

The matrices M ki and N k i have the same eigenvalues A. The 
advantage of the new method using Nh in (11) instead of M k i 
in (9) is that no differentiation of the operator D is required, so 
even numerically computed or tabulated values of u, a, and j8 
produce accurate results. 

2.4. Choice of bt and numerical handling of (11) 

For the set of base functions, we adopt the free magnetic de- 
cay modes whose magnetic fields ft, are known analytically in 
V + E in terms of the defining scalars P and T as described in 
Appendix A. The decay modes are continuous through dV and 
potential in E, so they satisfy the boundary conditions. They are 
characterised by three numbers, the radial order n, the latitudinal 
degree /, and the azimuthal order m. 

Another advantage of the decay modes is that they are self- 
adjoint on V + E so that the adjoint functions are the complex 
conjugates b k = b* k and likewise j k = j* k . Normalisation on V+E, 
i.e. (b* k ,bi)v+E = (j k > a i)v = 8 k h is thus readily achieved, see 
Appendix A. 

The computation of the matrix elements A^, is now straight- 
forward. Once we know the matrix elements, the eigenvalue 
problem (11) is solved numerically using LAPACK routines 
(http://www.netlib.org/lapack), and we obtain the eigenvalues A k 
and eigenvectors \e ki }, such that 

B k = e ki bi (12) 
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Table 1. Eigenvalues of the fundamental dipolar mode and the fifth and 
tenth overtones of the ^-sphere. 



"max 


A 


As 


'llO 


4 


-0.0241 


-116.936 




8 


-0.0034 


-115.110 


-349.489 


12 


-0.0010 


-115.058 


-349.097 


16 


-0.0004 


-115.047 


-349.050 


20 


-0.0002 


-115.043 


-349.036 


24 


-0.0001 


-115.041 


-349.031 


28 


-0.0001 


-115.040 


-349.028 


32 


-0.0000 


-115.040 


-349.026 


36 


-0.0000 


-115.039 


-349.026 


HvG93 


0.0000 


-115.04 


-349.02 



is eigenfunction ofVxD with eigenvalue A k . Each mode k con- 
tains, in general, a mixture of n, I, and m values. 

In the following we consider only velocities and dynamo co- 
efficients that are independent of azimuth ip, but this is not a 
necessary constraint. Thus each value of m can be treated sep- 
arately. Although we present only results for m — here, we 
have tested and applied other values of m as well. We employ 
the robust Gauss-Legendre quadrature in r and cos i? to compute 
the matrix elements since the basis functions are heavily oscilla- 
tory in r for high values of n and in 9 for high degree I. For the 
Gauss-Legendre integration we used 66 quadrature points here 
in the radial and 80 in the latitudinal direction, respectively. In 
general, this depends of course on the required resolution. 

2.5. Adjoint eigenfunctions 

We now show how one may construct the adjoint set of eigen- 
functions B p of a set of eigenfunctions B, of the dynamo opera- 
tor VxD. Although these adjoints are not needed in the present 
paper, they appear in applications. For example, let B be the 
actual magnetic field of the dynamo, then it is often advanta- 
geous to expand B in dynamo eigenfunctions, i.e., B(r,t) = 
Xj Ci(t)Bi(r). To find the coefficients c,(0 we use the adjoint set, 
to find c, = (B h B) v+E = (/,-, A) v . 

This illustrates that we need the adjoints B p , and these may 
be constructed as follows. Let B k = e ki bi be the representation of 
B k in terms of the self-adjoint magnetic decay modes as above. 
Then we write B p = f* pi b* and J p = f*J*, and we require 

S pk = (J p , A k ) v = if*J% e kj aj) v = f pi e k j5ij = f pi e ki . (13) 

Therefore f 1 = e _1 in matrix notation, and we find a unique 
biorthogonal set. Here, f indicates the Hermitean adjoint, ft = 
(f r )*, where T indicates the transposed and * complex conjuga- 
tion. 

Three important messages follow from this construction: (i) 
the adjoint of V x B is V x B, that is, the adjoint operation com- 
mutes with V; (ii) to obtain the adjoint eigenfunctions, it is not 
necessary to know the explicit form of the adjoint dynamo oper- 
ator VxD; and (iii) the eigenfunctions and their adjoints have the 
same boundary conditions because they are a linear combination 
of the decay modes and their complex conjugates, respectively. 



Table 2. Eigenvalues of the fundamental mode for the tiS 2 flow for two 
magnetic Reynolds numbers R m . 



"max — Anax 


(R m = 10) 


(Rm = 100) 


4 


-8.02911 


-5.88075 


8 


-8.02625 


-6.84495 


12 


-8.02626 


-6.92613 


16 


-8.02627 


-6.92859 


20 


-8.02627 


-6.92866 


24 


-8.02627 


-6.92869 


LJ05 


-8.01600 


-6.92885 



3. Test results 

3.1. a 2 -sphere 

We first consider the so-called ff 2 -sphere of unit radius r = 1, 
represented by u — B = 0, a t j = R a 6ij, and tj = 1, which also 
can be treated analytically (Krause & Radler 1980, Chap. 14). 
The eigenvalues are independent of azimuth m, and the eigen- 
functions decouple in latitudinal quantum number I. For R a + 0, 
the modes couple in radial number n, as they do between the 
poloidal and toroidal components. 

For R a = 4.493409458, 1 the first mode is a stationary dipole, 
while the overtones decay with the rates given by Hoyng & van 
Geffen (1993) (HvG93). We successfully reproduced the funda- 
mental mode and the overtones. In Table 1 we consider the con- 
vergence in the eigenvalues as a function of the maximum radial 
number n max for some dipolar (Z = 1) modes. Higher I modes 
behave similarly. We also reproduced the eigenfunction plots as 
provided by Krause & Radler (1980). 

3.2. Spherical flows 

As a next test, we apply the spherical stationary t\S2 (MDJ) flow 
of Livermore & Jackson (2004) which is given by 

u = u K- l Vx(r 2 (l -r 2 )P?(cos#)f) 

+ M o^'eVxVx(r 3 (l - r 2 ) 2 P° 2 (cos ff)rj (14) 

with K- { = V9009/572 and e = 0.5 V143/1008 such that the 
rms poloidal to toroidal energy ratio is 0.5, and the flow has an 
rms value of uq. P'" and f are defined in Appendix A. 

Like Livermore & Jackson (2005) we consider the axisym- 
metric {m = 0) and equatorially antisymmetric magnetic field so- 
lution for a unit sphere (ro = 1) embedded in a vacuum. Table 2 
shows the convergence in the eigenvalue with the largest real 
part as a function of truncation n max and / max for two magnetic 
Reynolds numbers R m = u Q r /T] =10 and R m = 100, together 
with the converged values given by Livermore & Jackson (2005) 
(LJ05). There is a difference of about one permille between their 
value and ours foxR m = 10. 

3.3. a 2 and a 2 D.-dynamos 

We reproduced the critical dynamo numbers R„ further for the 
dipolar (Z = 1) mode of an isotropic a 2 -dynamo with a rr = 
<*§§ = <X(p<p = R a sm Nn(r - r,) and N = 1 , 2 in a spherical 
shell of inner and outer radius r, and ro with ro - r, ■ — 1 and 



1 The numerical value of R a is equal to the first zero of the spherical 
Bessel function j\ . 
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Table 3. Eigenvalues of the first and fifth modes of the a 2 £l-dynamo (see Sect. 3.3). 



New method Classical method 



^max 


^max 


*o (S) 




-MA) 


M (S) 




M (A) 


8 


8 


(16.572,0.0) 


(- 


-25.767, ±45.538) 


(16.577,0.0) 


(- 


-25.763, ±45.529) 


12 


12 


(16.053,0.0) 


(- 


-33.126, ±39.444) 


(16.053,0.0) 


(- 


-33. 125, ±39.443) 


16 


16 


(16.052,0.0) 


(- 


-33.201, ±40. 168) 


(16.052,0.0) 


(- 


-33.201, ±40.168) 


20 


20 


(16.053,0.0) 


(- 


-33.200, ±40.157) 


(16.053,0.0) 


(- 


-33.200, ±40. 157) 


8 


20 


(16.052,0.0) 


(- 


-33.203, ±40. 158) 


(16.057,0.0) 


(- 


-33.197, ±40.163) 


12 


20 


(16.052,0.0) 


(- 


-33.200, ±40.157) 


(16.052,0.0) 


(- 


-33.199, ±40.157) 



r,/ro = 0.35 and 0.8 surrounded by a vacuum and either an in- 
sulating or a conducting inner core, as reported in Table 2 of 
Schubert & Zhang (2001). With this test we treated in particu- 
lar two different aspect ratios of a spherical shell (a thick and a 
thin one) and two different molecular diffusivities (insulating or 
conducting) of the inner core. 

Finally we applied our method to an a 2 f2-dynamo of Jiang & 
Wang (2006) who employ the classical eigenvalue treatment for 
the poloidal and toroidal scalars P and T expanded in spherical 
harmonics in the angular coordinates and in Chebychev polyno- 
mials in r-direction. We set u r — u# = 0, = Rnr 3 sin 3 ft, aij = 
Ptjk = 0, excepta^ = a§§ = a w = R a sin27r(r-r,)/(ro-r,) cos i? 
with r, = 0.5, ro = 1, embedded in a vacuum inside and outside, 
and R a = aoro/t] = 10 and R& = uoro/r] = 1000. Some re- 
sults obtained by the new and the classical methods are compiled 
in Table 3. Numbers in parentheses (. . . , . . .) are the real and 
imaginary parts of complex eigenvalues. The real part denotes 
the growth rate, the imaginary part the frequency of the mode 
in units of rj/rl 1 . The modes are axisymmetric (m = 0), the fun- 
damental mode is monotonously growing and symmetric (indi- 
cated by S ) with respect to the equator, and the fourth overtone is 
damped, oscillatory and antisymmetric (indicated by A). Modes 
with higher m are more strongly damped. n max refers to the max- 
imum radial number of the decay modes (spherical Bessel func- 
tions) for the new method and to the maximum degree of the 
Chebychev polynomials for the code of Jiang & Wang (2006), 
respectively. Since the modes have smaller length scales in lat- 
itudinal than in radial direction, higher values of I than of n are 
required for convergence. We find remarkably similar conver- 
gence of the eigenvalues for both methods. This also applies to 
modes with higher m. Of course we have also verified that the 
eigenfunctions obtained with the two methods are identical. 



4. Geodynamo models 

Having proven that the new method works correctly and ef- 
ficiently, we now apply it to determine the eigensolutions 
of the dynamo operator with mean-field coefficients obtained 
from self-consistent numerical simulations of the geodynamo. 
For a recent review of numerical geodynamo simulations, see 
Christensen & Wicht (2007). Schrinner et al. (2007) developed 
an efficient method of calculating all tensorial mean-field coef- 
ficients a and fi and compared the results of mean-field and di- 
rect numerical simulations of the geodynamo. We plan to use the 
eigenmodes of the dynamo equation to decompose the magnetic 
field of the numerical simulations and to determine the statistical 
properties of the mode coefficients (Hoyng 2009) to analyse the 
working of the geodynamo. 



4.1. Benchmark dynamo 

We examine a quasi-steady geodynamo model which has been 
used before as a numerical benchmark dynamo (Christensen et 
al. 2001, case 1). The governing parameters are Ekman num- 
ber E = 10~ 3 , Rayleigh number Ra = 100, Prandtl number 
Pr — 1, and magnetic Prandtl number Pm = 5. The convection 
pattern is columnar with a natural 4-fold azimuthal symmetry 
and is stationary except for an azimuthal drift. The intensity of 
the fluid motion is characterised by a magnetic Reynolds num- 
ber of R m =! 40, defined with a characteristic flow velocity, the 
thickness of the convecting shell, and the molecular magnetic 
diffusivity. The magnetic energy density exceeds the kinetic one 
by a factor of 20. 

In Schrinner et al. (2007), the mean-field coefficients are 
derived from the numerical simulation. We solved the dynamo 
equation with these mean-field coefficients by the new method 
and obtained the eigenvalues and eigenfunctions. Since the co- 
efficients are spatially variable to a considerable degree, con- 
verged solutions require high truncation levels in n and /. The 
eigenvalues of the first two modes are shown in Table 4. Beyond 
"max - 20 and Z max 20, the eigensolution of the first mode does 
not change significantly and is displayed in Fig. 1. The conver- 
gence of the second mode requires a larger / max of about 32. The 
results for high values of n max may be affected by the spatial vari- 
ation of the mean-field coefficients and would require more than 
66 radial quadrature points to compute the matrix elements. 

A comparison of Fig. 1 with its counterpart Fig. 10 of 
Schrinner et al. (2007) shows that the field of the antisymmetric 
fundamental mode resembles the field of an initial-value mean- 
field dynamo calculation remarkably well as it does the axisym- 
metric component of the direct numerical simulation. The mode 
here grows slightly with a rate around Aq 4.2^/L 2 , the field 
of the initial value calculation decays slightly with a rate of ap- 
proximately -0.25r]/L 2 , while the solution of the direct numeri- 
cal simulation is stationary 2 . Here L = ro - r, = 1 is the thickness 
of the spherical shell. The difference in these rates between the 
eigenvalue and initial value calculation comes from the higher 
numerical diffusivity of the latter at the chosen resolution of 33 
radial and 80 latitudinal grid points. The difference is actually 
small, much less than one effective decay rate, because the rel- 
evant turbulent diffusivity, described by the /3 coefficient with 
values up to 33^, is much higher than the molecular one. 

Besides the true physical eigenmodes, we find growing 
unphysical spurious eigenmodes. Their eigenvalues depend 
strongly on the resolution, and their eigenfunctions are highly 
structured. We attribute their appearance to a locally confined in- 
appropriate parametrisation of the mean electromotive force by 



In Schrinner et al. (2007) the mean flow entered with a sign er- 
ror into the initial value calculation, leading to a stronger decay of 
-3.5rj/L 2 . We apologise and correct this value here. 
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Table 4. Eigenvalues in units of rj/L 2 of the first two eigenmodes of the 
benchmark dynamo. 



'^max 


^max 


A (A) 


A l (5) 


12 


12 


+4.960 


-8.605 


16 


16 


+4.235 


-8.110 


20 


20 


+4.180 


-8.362 


24 


24 


+4.195 


-6.620 


28 


28 


+4.255 


-7.275 


32 


32 


+4.275 


-6.015 


16 


32 


+4.382 


-6.777 



Notes. The symmetry with respect to the equator is marked A for an- 
tisymmetric and S for symmetric. The eigensolution marked bold is 
displayed in Fig. 1. 




Fig. 1. Magnetic field structure of the fundamental antisymmetric eigen- 
mode for the benchmark dynamo. Compare with Fig. 10 of Schrinner et 
al. (2007). For each plot the grey scale is separately adjusted to its maxi- 
mum modulus with white as negative and black as positive. The contour 
lines correspond to ±0.1, ±0.3, ±0.5, ±0.7, and ±0.9 of the maximum 
modulus. 



the mean-field coefficients a and J3 (Schrinner et al. 2007). The 
spurious modes are present neither in the initial value calcula- 
tion nor in the following example of a time-dependent dynamo, 
because of a higher numerical and molecular diffusivity, respec- 
tively. 

4.2. A time-dependent dynamo in the columnar regime 

The next example has stronger forcing with parameters E = 
10~ 4 , Ra = 334, Pr = 1, and Pm = 2. The numerical simulation 
by Olson et al. (1999, case 2) shows a highly time-dependent, but 
still dominantly columnar convection characterised by a mag- 
netic Reynolds number of R„, ^ 88. The magnetic energy ex- 
ceeds the kinetic energy by a factor of three. The magnetic 
field has a strong axial dipole contribution. Although chaotically 
time-dependent, the velocity field is symmetric and the magnetic 
field antisymmetric with respect to the equatorial plane. 

The mean-field coefficients are obtained as before by the 
test-field method of Schrinner et al. (2007). The coefficients are 
now of course also highly time-dependent. A time average yields 
coefficients that roughly resemble those for the benchmark dy- 
namo, although there are differences in some profiles and ampli- 
tudes. 



Table 5. Eigenvalues in units of rj/L 2 of the first two antisymmetric 
eigenmodes of the temporally averaged dynamo operator obtained from 
the time-dependent dynamo (case 2, Sect. 4.2). 



^raax 


^max 




A 




^3 


12 


12 


(- 


-4.520,0.0) 


(- 


35.046, ±10. 118) 


16 


16 


(- 


•4.278 , 0.0) 


(- 


35.090, ±10.256) 


20 


20 


(- 


-4.112,0.0) 


(- 


34.578, ±10.096) 


24 


24 


(- 


-3.930,0.0) 


( 


-34.770, ±9.874) 


28 


28 


(-3.880,0.0) 


(- 


34.766, ±10.236) 


32 


32 


(- 


-3.868,0.0) 


(- 


34.804, ±10.318) 


16 


32 


(- 


-3.874,0.0) 


(- 


34.830, ±10.310) 



Notes. The eigensolutions marked bold are displayed in Fig. 2. 



B,: A=(-3.87,0.00) Re(BJ: A=(-34.83, 10.31) lm(B,): A=(-34.83,10.31) 

Max: 1.44 Max: 0.58 Max: 1.37 




Fig. 2. Radial component of the first two antisymmetric eigenmodes for 
the case2 dynamo. Left: fundamental mode; middle: real part of the 
first overtone; right: imaginary part of the first overtone. Grey scales 
and contours as in Fig. 1. 



For the time-averaged dynamo operator the eigenvalues of 
the first two antisymmetric eigenmodes for various values of 
«max and Z max are shown in Table 5. It seems that a value of « max ^ 
16 is sufficient for convergence, while l max ^ 32 is needed. 
Figure 2 shows the eigenfunctions of these modes. The eigen- 
solutions for Ai = (-6.298 , 0.0) and A 2 = (-28.712 , ±5.364), 
values for « max = 16 and Z max = 32, are symmetric with respect 
to the equator. 

An initial-value, mean-field dynamo calculation with the 
same mean velocity and dynamo coefficients shows a slighly 
decaying solution with a decay rate of approximately -5.9i]/L 2 
which is to be compared with the eigenvalue Aq ^ -3.87?7/L 2 of 
the fundamental mode. Again, the turbulent diffusivity exceeds 
the molecular one by a factor of up to 23 in this case. The differ- 
ence in the decay rates is therefore much less than one effective 
decay rate. As for the benchmark dynamo, the profile of the anti- 
symmetric fundamental mode is again remarkably similar to the 
solution of the initial value calculation and to the axisymmetric 
component of the direct numerical simulation. 

A decomposition of the actual magnetic field of the simula- 
tion by Olson et al. (1999, case 2) in eigenfunctions of the time- 
averaged dynamo operator, i.e., B(r,t) = £j Ci{i)Bi{r), shows 
that the antisymmetric fundamental mode contributes to about 
75 percent and, together with the first antisymmetric overtone 
(see Table 5 and Fig. 2), to about 85 percent of the total mag- 
netic energy. The variability in time of the magnetic field of the 
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direct numerical simulation is reflected in the variability of the 
expansion coefficients. More details are presented in Schrinner 
et al. (2009). 

5. Conclusions and outlook 

We presented a new method for computing the eigenvalues and 
eigenfunctions of the induction and the dynamo equation. The 
method is based on the biorthogonality of the adjoint electric 
current and the vector potential with an inner product defined 
by a volume integral over the fluid domain. The advantage of 
the method is that the velocity and dynamo coefficients do not 
have to be differentiated. The method is therefore well-suited for 
spatially strongly variable dynamo coefficients. 

We tested the new method against the classical treatment and 
proved that it works correctly and efficiently. We applied it to 
two cases with dynamo coefficients derived from direct numeri- 
cal simulations of the geodynamo. The obtained dynamo eigen- 
modes are promising candidates for decomposing the magnetic 
field of the numerical simulations and for analysing the statis- 
tical properties of the mode coefficients as proposed by Hoyng 
(2009). 
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Appendix A: Free magnetic decay modes in a 
sphere or spherical shell embedded in vacuum 

We decompose the magnetic field in its poloidal and toroidal 
components 



B = VxVxPr + Vx7r 



(A.l) 



with defining scalars P(r, §, <p) and T(r, i?, <p) and unit vector in 
radial direction f = (1,0, 0) in spherical coordinates (r,ft,<p). 
The equation of free magnetic decay 



dB 

— = -nVxVxB 

at 

with constant magnetic diffusivity 77 then reads 

^) = ^ + -)(P,r) 

where A H is the horizontal Laplacian 

1 d . „ d 1 d 2 

■ sin 17 — + 



(A.2) 



(A.3) 



(A.4) 



r 2 sin 2 #d& d» r 2 sin 2 dip 2 ' 
The solutions are the free magnetic decay modes 

P lmn = exp (Alt) xMp ln x) 17(0, if) , (A.5) 
T lmn = exp(Alt)x gl (ti„x)Y™(V,<p) (A.6) 

with x = r/ro where ro is the radius of the sphere. The growth 
rates are given by 



< = -VPL/iq and An = -V&/>o 



(A.7) 



and are independent of the azimuthal degree m. The constants 
pi„ and ti n are 



Pin = jl-l,n and ti n = jl,n 



(A.8) 



where ji „ is the n-th zero of j/. The T™ are the spherical harmon- 
ics and normalised to unity by taking 



4n (l + m)\ 
21+1 (l-m)\ 



■1/2 



P™{cos&)e imf 



(A.9) 



using Ferrer's definition of the Legendre functions of first kind 
P"' with degree I and order m. 

For a sphere embedded in vacuum the radial functions are 
given by 



, ai„ jl(pinX) < X < 1 

MPl " X) - < a ln Hp ln )x- 1 x>\ 



gi(ti„x) 



bin jlitlnX) < X < 1 







X > 1 



(A. 10) 
(All) 
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with the spherical Bessel functions of first kind j/. This ensures 
regularity in the origin of the sphere, vanishing toroidal compo- 
nent at its outer boundary and smooth transition of the poloidal 
component to a potential field in the vacuum outside. 

For a spherical shell with inner radius r, (x, = r,-/ro) and 
outer radius ro (xq = 1) embedded in vacuum the radial functions 
inside the shell are given by 

fl(PlnX) = jl(pinX) ~ yi(PlnX)jl + l(pl„Xi)/y M (pi„Xi) (A. 12) 

and 

gl(tlnX) = jl(tl„x) - yi(tl n X)ji(ti n Xi)/yi(ti n Xi) , (A. 13) 

and the constants in the arguments are the roots of 

jl+\{PlnXi)yi-\(pi„) - jl-l(pin)yi+\(PlnXi) = (A.14) 

for pi„ and of 

jl(tln)yi(tl„Xi) - ji(ti n x,)yi(t ln ) = (A. 15) 

for ti n . Here yi are the spherical Bessel functions of second kind. 

The magnetic field of the decay modes B, is obtained by in- 
serting the spatial parts of the defining scalars P/ mn and T[ mn , re- 
spectively, into (A.l). Here we have comprised the three indices 
into one. The decay modes are self-adjoint on V + E, so that the 
adjoint functions are obtained simply by complex conjugation: 
Bk = B* k and likewise = J*,. Normalisation on V + E, i.e., 
(B* k , Bj) v +E = ( J* k , Ai) v = 8ku is thus straightforward. For a unit 
sphere the radial functions are normalised to unity by scaling the 
// with 

am = \^rKl + l)jl Ul fi(ji-un)} (A.16) 
and the gi with 

frz» = (y/(/+l)7ti(A«)) ^ • (A. 17) 

For a spherical shell the normalisation constants are more 
lengthy expressions, which we suppress here. 

The free magnetic decay modes form a complete and orthog- 
onal set of functions, and they obey the boundary conditions of 
the magnetic field between the dynamo volume V and the exte- 
rior vacuum E. 

We mention for completeness that the poloidal decay modes 
are not self-adjoint on V, i.e., (B* k , B,)v + 6u- If we like to work 
with an inner product defined on V, the adjoint functions B^ can 
be constructed by requiring (Bk,Bi)v = 5ki, similar to the one 
described in Sect. 2.5. 
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